home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / mrg.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  4.0 KB  |  150 lines
  1. /* rng/mrg.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_rng.h>
  23.  
  24. /* This is a fifth-order multiple recursive generator. The sequence is,
  25.  
  26.    x_n = (a_1 x_{n-1} + a_5 x_{n-5}) mod m
  27.  
  28.    with a_1 = 107374182, a_2 = a_3 = a_4 = 0, a_5 = 104480 and m = 2^31-1.
  29.  
  30.    We initialize the generator with x_n = s_n MOD m for n = 1..5,
  31.    where s_n = (69069 * s_{n-1}) mod 2^32, and s_0 = s is the
  32.    user-supplied seed.
  33.  
  34.    NOTE: According to the paper the seeds must lie in the range [0,
  35.    2^31 - 2] with at least one non-zero value -- our seeding procedure
  36.    satisfies these constraints.
  37.  
  38.    We then use 6 iterations of the generator to "warm up" the internal
  39.    state.
  40.  
  41.    With this initialization procedure the theoretical value of
  42.    z_{10006} is 2064828650 for s = 1. The subscript 10006 means (1)
  43.    seed the generator with s = 1, (2) do the 6 warm-up iterations
  44.    that are part of the seeding process, (3) then do 10000 actual
  45.    iterations.
  46.  
  47.    The period of this generator is about 2^155.
  48.  
  49.    From: P. L'Ecuyer, F. Blouin, and R. Coutre, "A search for good
  50.    multiple recursive random number generators", ACM Transactions on
  51.    Modeling and Computer Simulation 3, 87-98 (1993). */
  52.  
  53. static inline unsigned long int mrg_get (void *vstate);
  54. static double mrg_get_double (void *vstate);
  55. static void mrg_set (void *state, unsigned long int s);
  56.  
  57. static const long int m = 2147483647;
  58. static const long int a1 = 107374182, q1 = 20, r1 = 7;
  59. static const long int a5 = 104480, q5 = 20554, r5 = 1727;
  60.  
  61. typedef struct
  62.   {
  63.     long int x1, x2, x3, x4, x5;
  64.   }
  65. mrg_state_t;
  66.  
  67. static inline unsigned long int
  68. mrg_get (void *vstate)
  69. {
  70.   mrg_state_t *state = (mrg_state_t *) vstate;
  71.  
  72.   long int p1, h1, p5, h5;
  73.  
  74.   h5 = state->x5 / q5;
  75.   p5 = a5 * (state->x5 - h5 * q5) - h5 * r5;
  76.   if (p5 > 0)
  77.     p5 -= m;
  78.  
  79.   h1 = state->x1 / q1;
  80.   p1 = a1 * (state->x1 - h1 * q1) - h1 * r1;
  81.   if (p1 < 0)
  82.     p1 += m;
  83.  
  84.   state->x5 = state->x4;
  85.   state->x4 = state->x3;
  86.   state->x3 = state->x2;
  87.   state->x2 = state->x1;
  88.  
  89.   state->x1 = p1 + p5;
  90.  
  91.   if (state->x1 < 0)
  92.     state->x1 += m;
  93.  
  94.   return state->x1;
  95. }
  96.  
  97. static double
  98. mrg_get_double (void *vstate)
  99. {
  100.   return mrg_get (vstate) / 2147483647.0 ;
  101. }
  102.  
  103.  
  104. static void
  105. mrg_set (void *vstate, unsigned long int s)
  106. {
  107.   /* An entirely adhoc way of seeding! This does **not** come from
  108.      L'Ecuyer et al */
  109.  
  110.   mrg_state_t *state = (mrg_state_t *) vstate;
  111.  
  112.   if (s == 0)
  113.     s = 1;    /* default seed is 1 */
  114.  
  115. #define LCG(n) ((69069 * n) & 0xffffffffUL)
  116.   s = LCG (s);
  117.   state->x1 = s % m;
  118.   s = LCG (s);
  119.   state->x2 = s % m;
  120.   s = LCG (s);
  121.   state->x3 = s % m;
  122.   s = LCG (s);
  123.   state->x4 = s % m;
  124.   s = LCG (s);
  125.   state->x5 = s % m;
  126.  
  127.   /* "warm it up" with at least 5 calls to go through
  128.      all the x values */
  129.  
  130.   mrg_get (state);
  131.   mrg_get (state);
  132.   mrg_get (state);
  133.   mrg_get (state);
  134.   mrg_get (state);
  135.   mrg_get (state);
  136.  
  137.   return;
  138. }
  139.  
  140. static const gsl_rng_type mrg_type =
  141. {"mrg",                /* name */
  142.  2147483646,            /* RAND_MAX */
  143.  0,                             /* RAND_MIN */
  144.  sizeof (mrg_state_t),
  145.  &mrg_set,
  146.  &mrg_get,
  147.  &mrg_get_double};
  148.  
  149. const gsl_rng_type *gsl_rng_mrg = &mrg_type;
  150.